Medford Disaster Router¶
Michael Ruggiero
In [1]:
# Function imports
import numpy as np
import pandas as pd
import networkx as nx
import osmnx as ox
import geopandas as gpd
import matplotlib.pyplot as plt
import matplotlib.cm as cm
import matplotlib.colors as colors
ox.config(use_cache=True, log_console=True)
ox.__version__
%config InlineBackend.figure_format = 'retina'
%matplotlib inline
1. Problem Statement ¶
The recent Merrimack Valley gas explosions reveled the need for cities to quickly deploy first responders while accounting for damaged roads. However, information during emergencies can be incomplete and properly dispatching police requires multiple information streams.
This project attempts to address this problem with the following workflow:
(1) Incorporate DOT api information about construction and accidents along flagged roads.
(2) Incorporate Twitter and Google information to estimate traffic patterns.
(3) Simulate a disaster zone in Medford that destroys a percentage of roads and creates emergency locations.
(4) Optimize dispatch routes for police to reach each disaster location.
2. Data Collection/Import Medford¶
Basic information about an importated city. Due to my site specific knowledge, I decided to protoype this model on the city of Medford. Other city disasters will be investigated. I-93 was removed from the road network to simplify the initial analysis.¶
In [2]:
# get a graph for some city
place = {'city' : 'Medford',
'state' : 'MA',
'country' : 'USA'}
medford = ox.graph_from_place(place, network_type='bike')
fig, ax = ox.plot_graph(medford)
In [14]:
# what sized area does our network cover in square meters?
medford_proj = ox.project_graph(medford)
nodes_proj = ox.graph_to_gdfs(medford_proj, edges=False)
graph_area_m = nodes_proj.unary_union.convex_hull.area
print(graph_area_m)
# show basic stats about the network
medford_stats = ox.basic_stats(medford_proj,
area=graph_area_m,
clean_intersects=True,
circuity_dist='euclidean')
# medford_stats = ox.extended_stats(medford_proj, ecc=True, bc=True, cc=True)
pd.DataFrame(medford_stats).set_index("streets_per_node_counts")
Out[14]:
In [15]:
# Edge and Node projection
nodes_med, edges_med = ox.graph_to_gdfs(medford, nodes=True, edges=True)
In [16]:
nodes_med.head(2)
Out[16]:
In [17]:
nodes_med.shape
Out[17]:
2a. Major Roads Functions¶
To avoid cluttering up the notebook, all major functions have been exported into the "graph_functions.py" python file. See this file for futher information.¶
In [18]:
#Functions built for this project
# %load graph_functions.py
import graph_functions as gf
The "node_roader" function returns a data dictionary of the essential information from a plotted object.¶
In [19]:
medford_data = gf.node_roader(medford)
In [20]:
medford_data.keys()
Out[20]:
In [21]:
edges_med, nodes_med, medford = medford_data["edges"], medford_data["nodes"], medford_data["graph"]
In [22]:
nodes_med.major_inter.value_counts()
Out[22]:
In [23]:
#Long/Lat information
nodes_med[["x","y"]].describe()
Out[23]:
In [25]:
nodes_med.head(2)
Out[25]:
2b. Routing¶
The "random point" and "get_nearest_node" functions are designed to pair any longitude and latitude location with a node on the road network. Routing is handled by the Networkx package.¶
In [26]:
#A function to generate a random lat/long point
gf.random_point()
Out[26]:
In [27]:
random_point_A = gf.random_point()
random_point_B = gf.random_point()
# get the nearest network node to each random point
orig_node = ox.get_nearest_node(medford, random_point_A)
dest_node = ox.get_nearest_node(medford, random_point_B)
In [28]:
# find the route between these nodes then plot it
route = nx.shortest_path(medford, orig_node, dest_node, weight='length')
fig, ax = ox.plot_graph_route(medford, route, node_size=0)
In [29]:
# Length of route in meters
nx.shortest_path_length(medford, orig_node, dest_node, weight='length')
Out[29]:
In [30]:
# Absolute Distance between two nodes (as the crow flies)
ox.great_circle_vec(medford.node[orig_node]['y'], medford.node[orig_node]['x'],
medford.node[dest_node]['y'], medford.node[dest_node]['x'])
Out[30]:
2c. Disaster Radius¶
"Disaster_generator" creates a graph object around a node given a radius. This subgraph serves as the region effected by the disaster.¶
In [31]:
d_location = nodes_med[["y","x"]].sample(1).values[0]
In [33]:
print(d_location)
In [32]:
#returns a graphs object of a disater
disaster = gf.disaster_generator(nodes_med, location_point = d_location, radius = 2500)
In [34]:
nodes_dis, edges_dis = ox.graph_to_gdfs(disaster)
The disaster is then reimposed onto the original city graph.¶
In [35]:
#Set Colors
ec = ['lightcoral' if i in disaster.edges() else 'green' for i in medford.edges()]
nc = ['red' if i in disaster.nodes() else 'blue' for i in medford.nodes()]
#Plot energency grid
ox.plot_graph(medford, node_size=15, node_color = nc ,edge_color=ec)
Out[35]:
In [36]:
nx.degree_histogram(disaster)
Out[36]:
The "road_kill" function destroys a given percentage of roads inside the disaster zone. These destroyed roads are removed from edges of the city network. Comparing the images above and below will demonstrate that some nodes are now completely isolated as all connecting roads are destroyed.¶
In [41]:
medford_disaster = gf.road_kill(disaster,
nodes_med,
edges_med,
kill_percentage = .8 )
2d. City Visualizations¶
Various city road information and city information¶
In [ ]:
# get the boundary polygons for neighboring cities, save as shapefile, project to UTM, and plot
place_names = ['Medford, MA, USA',
'Arlington, MA, USA',
'Somerville, MA, USA',
'Malden, MA, USA',
# 'Melrose, MA, USA',
'Winchester, MA, USA',
'Stoneham, MA, USA',
'Everett, MA, USA',
]
medford_area = ox.gdf_from_places(place_names)
fig, ax = ox.plot_shape(medford_area)
In [ ]:
# highlight one-way roads
oneway = ox.graph_from_place(place, network_type='drive')
ec = ['r' if data['oneway'] else 'b' for u, v, key, data in oneway.edges(keys=True, data=True)]
fig, ax = ox.plot_graph(oneway, node_size=0, edge_color=ec, edge_linewidth=1.5, edge_alpha=0.5)
In [ ]:
#Neighboor hoods buffered
neigboors_buffered = ox.gdf_from_places(place_names, gdf_name='neighboors', buffer_dist=250)
fig, ax = ox.plot_shape(neigboors_buffered, alpha=0.7)
3. Data Creation¶
"Random_zone_picker" creates random crosstown traffic .CSVs for google to analyze.¶
In [42]:
# DF builder of start/end nodes with long/lat and times
# Already done, commenting out to avoid clutter in notebook
router = gf.random_zone_picker(nodes_med, 10)
#router.to_csv("more_zones.csv")
In [44]:
router.head(3)
Out[44]:
Police are not restricted to one city, thus, a region around Medford is mapped to potentially use as routes around the disaster. Here a 5.5 km radius is plotted around city hall.¶
In [45]:
# network from address, including only nodes within 5.5km along the network from city hall
neighbor = ox.graph_from_address(address='85 George P. Hassett Drive,\
Medford, MA 02155',
distance=5500,
distance_type='network',
network_type='walk')
# you can project the network to UTM (zone calculated automatically)
neighbor_projected = ox.project_graph(neighbor)
In [46]:
# Edge and Node projection
nodes_area, edges_area = ox.graph_to_gdfs(neighbor_projected,
nodes=True,
edges=True)
nc = ['blue' if i in medford.nodes() else 'green' for i in neighbor.nodes()]
ox.plot_graph(neighbor,
node_size=15,
node_color = nc);
In [47]:
neighbor_data = gf.node_roader(neighbor)
In the above coloration, black nodes serve as significant intersection between major roads.¶
In [50]:
black = neighbor_data["nodes"]
black[black.color == "black"].head(2)
Out[50]:
In [51]:
neighbor_data.keys()
Out[51]:
In [52]:
nodes_area, edges_area, neighborhood = neighbor_data["nodes"], neighbor_data["edges"], neighbor_data["graph"]
Major node connection vector¶
In [54]:
for i in neighbor_data['major_intersections']:
print("Number of major nodes: ",
len(neighbor_data['major_intersections'][i]),
"\t with degree ",
i,)
In [55]:
neighbor_data['major_intersections'].keys()
Out[55]:
In [56]:
# find the route between these nodes then plot it
route = nx.shortest_path(neighbor, 3809869827, 65339393, weight='length')
fig, ax = ox.plot_graph_route(neighbor, route, node_size=0)
Checks to see if a neighborhood node is the city and creates a column based on data in the node geopandas dataframe¶
In [57]:
_ = gf.city_node(nodes_med,"medford",nodes_area)
In [58]:
nodes_area.head()
Out[58]:
4 Google Directions Import¶
The following chart shows the approximation between the intersection use and Pareto Principle. 80% of the traffic is concentrated in 20% of the intersections.¶
In [59]:
gmap, complete_node, node_occurances = gf.directionsToDataframe("10.Google_directions_api_all6k.csv")
4a. Google top 10 nodes¶
In [61]:
gmap.head()
Out[61]:
In [64]:
google_top = gf.traffic_node(nodes_area, node_occurances)
g_list = list(google_top[google_top.traffic_importance > 887].osmid.values)
In [65]:
google_top[google_top.traffic_importance > 887]
Out[65]:
The following graph shows the location of the 10 most important nodes for traffic. From personal experience, this information is correct.¶
In [66]:
# Edge and Node projection
nodes_area, edges_area = ox.graph_to_gdfs(neighbor_projected,
nodes=True,
edges=True)
g_size = [25 if i in g_list else 0 for i in neighbor.nodes()]
ox.plot_graph(neighbor,
node_size=g_size,
node_color = "red", edge_alpha= .5);
5. Initial Monte Carlo Simulation¶
Rules:¶
- Simulation
- All disasters originate in Medford with no spill over
- 5 patrol officers who maybe inside disaster radius
- 5 officers dispatched from station
- Police do not follow one way road rules
- 10 randomly selected node emergencies inside disaster radius
- Best route chosen for all officers
- Officers dispatched for on foot for unreachable emergencies (distance 10,000 listed)
- 10,000 simulations
- I-93 is taken out by the disaster
- Disaster radius randomly selected
- Questions to answer
- Average police travel distance
- Most used nodes by police
In [68]:
import time
import pickle
from scipy.optimize import linear_sum_assignment
#Populate police station
"""
Medford Police:
100 Main St, Medford, MA 02155
42.415811 -71.110152
"""
police_node = gf.nodefinder(nodes_med,
lat = 42.415811,
lng = -71.110152)
def monte_simulator(p_e_samples, police_location, simulations, name = ""):
main_start = time.time()
police_node = gf.nodefinder(nodes_med,
lat = police_location[0],
lng = police_location[1])
police_station = [police_node] * p_e_samples
p_e_samples *= 2
monte_simulation = {}
#Iteration though simulation
for sim in range(0,simulations):
#Initialization
start = time.time()
print("Starting simulation: {} at time: {}".format(sim,round(start - main_start,2)))
monte_dis = {}
#Generates disaster and radius
disaster_location = nodes_med.sample(1)
monte_dis["dis_location"] = disaster_location
monte_dis["dis_rad"] = np.random.uniform(low = .4, high = .9) * 4000
monte_dis["disaster"] = gf.disaster_generator(nodes_med,
location_point = disaster_location[["y","x"]].values[0],
radius = monte_dis["dis_rad"],
plotter = 0)
monte_dis["nodes_effected"] = len(monte_dis["disaster"].nodes())
monte_dis["roads_effected"] = len(monte_dis["disaster"].edges())
#Randomly destroys roads in disaster radius
monte_dis["remain_per"] = np.random.uniform(low = .35, high = 1)
monte_dis["remaining"] = gf.road_kill(monte_dis["disaster"],
nodes_area, #These come from the neighborhood
edges_area,
kill_percentage = monte_dis["remain_per"],
plotter = 0)
#Make dictionaries for emergencies and police
e_list = np.random.choice(list(monte_dis["disaster"].nodes()),
replace=False,
size = p_e_samples)
p_list = list(nodes_med.sample(int(p_e_samples / 2)).osmid.values)
p_list.extend(police_station)
monte_dis["emergency"] = {}
monte_dis["patrol"] = {}
for i in range(0,p_e_samples):
monte_dis["patrol"]["officer_{}".format(i)] = p_list[i]
monte_dis["emergency"]["emergency_{}".format(i)] = e_list[i]
rows = list(monte_dis["patrol"].keys())
columns = list(monte_dis["emergency"].keys())
#Distance Matrix
d_matrix = pd.DataFrame(columns = columns,
index = rows)
for i in range(len(rows)):
for j in range(len(columns)):
try:
path = nx.shortest_path_length(monte_dis["remaining"],
monte_dis["patrol"][rows[i]],
monte_dis["emergency"][columns[j]],
weight='length')
except:
path = 10_000
d_matrix.loc[rows[i], columns[j]] = path
monte_dis["distance_matrix"] = d_matrix.astype(float)
#Use Hungarian algorithm for linear sum assignment
lsa = linear_sum_assignment(monte_dis["distance_matrix"].to_numpy())
adj_matrix = np.matrix(np.zeros((10, 10)))
#optimized distances
od_list = []
for i in range(len(lsa[0])):
od = monte_dis["distance_matrix"].iloc[lsa[0][i],lsa[1][i]]
adj_matrix[lsa[0][i],lsa[1][i]] = od
od_list.append(od)
dispatch = {"officer_{}".format(lsa[0][i]):["emergency_{}".format(lsa[1][i]),
od_list[i]] for i in range(len(lsa[0]))}
monte_dis["ad_matrix"] = adj_matrix
monte_dis["dispatch"] = dispatch
monte_dis["dispatch_distance"] = sum(od_list)
monte_simulation[sim] = monte_dis
if sim % 200 == 0:
print("pickle dump")
pickle.dump(monte_simulation, open( "monte_simulation_{}_{}.p".format(name,sim), "wb" ))
end = time.time()
print("Ending simulation: {}. Total {} seconds for run\n".format(sim,round(end - start,2)))
return monte_simulation
Demo of five simulated disasters¶
In [69]:
monte_simulation = monte_simulator(5,[42.415811,-71.110152], 5, "demo")
Simulations inside data dictionary¶
In [70]:
monte_simulation.keys()
Out[70]:
Terms inside of each simulation¶
In [79]:
monte_simulation[0].keys()
Out[79]:
Two sample optimized dispatches with corresponding distance in meters¶
In [80]:
monte_simulation[0]['dispatch']
Out[80]:
In [82]:
monte_simulation[2]['dispatch']
Out[82]:
In [86]:
monte_simulation[2]["ad_matrix"]
Out[86]:
Data frame generation¶
In [71]:
size_list = range(len(monte_simulation))
In [73]:
df = pd.concat([monte_simulation[i]["dis_location"].copy() for i in size_list]).reset_index(drop=True)
In [76]:
extra_data = ['dis_rad', 'nodes_effected', 'roads_effected', 'remain_per', 'dispatch_distance']
In [77]:
for i in range(len(monte_simulation)):
for j in extra_data:
df.loc[i, j] = monte_simulation[i][j]
In [78]:
df.columns
Out[78]:
In [84]:
X = df[['osmid','x', 'y','dis_rad', 'nodes_effected', 'roads_effected',
'remain_per', 'dispatch_distance']]
In [85]:
X.head()
Out[85]:
In [ ]: